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Abstract. We investigate the late-time asymptotic behavior of solutions to 
nonlinear hyperbolic systems of conservation laws containing stiff relaxation 
terms. First, we introduce a Chapman-Enskog-type asymptotic expansion and 
derive an effective system of equations describing the late-time/stiff relaxation 
singular limit. The structure of this new system is discussed and the role 
of a mathematical entropy is emphasized. Second, we propose a new finite 
volume discretization which, in late-time asymptotics, allows us to recover a 
discrete version of the same effective asymptotic system. This is achieved pro- 
vided we suitably discretize the relaxation term in a way that depends on a 
matrix-valued free-parameter, chosen so that the desired asymptotic behavior 
is obtained. Our results are illustrated with several models of interest in con- 
tinuum physics, and numerical experiments demonstrate the relevance of the 
proposed theory and numerical strategy. 



1. Introduction 

Motivations. We are interested in the numerical approximation of solutions to 
nonlinear hyperbolic systems of conservation laws containing stiff relaxation terms. 
An extensive literature is available on such systems since they arise in many phys- 
ical problems of interest, for instance, in the modeling of complex multiphase 
flows involving phase transitions or kinetic-type phenomena. (See, for instance, 
[HI d31 HH dH 130] ■) Stiff relaxation source terms are essential to model phenomena 
involving distinct physical time-scales. The derivation and the analysis (existence, 
stability) of an effective system of equations (also of hyperbolic type as the origi- 
nal system) as the relaxation time goes to zero is required for understanding the 
behavior of general solutions. (See, for instance, [25l I3T1 l34l [35] .) 

In the present paper, we go beyond these classical works and investigate the 
late-time behavior of solutions to systems with stiff relaxation and, specifically, 
we consider the class of systems (e > 0) 

(1.1) edtU + d x F(U) = -^Q, t > 0, x e R, 

e 

where the main unknown U : M x M + — > ft takes its values in a convex set ft C 
Mr. The associated homogeneous first-order system — obtained by neglecting the 
relaxation term R : ft — > M. N in the right-hand side of (|l.ll) — is assumed to be 
hyperbolic in the sense that the N x N matrix A := DjjF admits real eigenvalues 
and a full basis of eigenvectors. 
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In comparison with classical works on this subject [TT] as well as also [U El [TQl 
[25] . the main novelty in the present work lies in the fact that the term edtU is 
rescaled to be proportional to e. When e tends to zero, this will lead us to limit 
solutions whose behavior is quite distinct from the ones of standard relaxation 
limits. Indeed, as we establish below, the effective problem associated with (II. ip 
turns out to be a diffusion problem, in which the diffusion operator is determined 
by the (nonlinear) relaxation term R, in a way explicitly determined in the present 
worlQ. The relaxation map R : fl — > R N is assumed to be sufficiently regular 
and satisfy the conditions introduced by Chen, Levermore, and Liu for the 
"hyperbolic to hyperbolic" relaxation problem. 

First of all, we assume the existence of a (constant) nx N matrix Q with (max- 
imal) rank n < N such that 

(1.2) QR{U) = 0, u en. 

From this, it follows that, if U is a solution to (fTl]), then QU : I x 1+ -> R N 
satisfies the system of n conservation laws 

(1.3) ed t (QU) + d x (QF(U)) = 0, 
in which the unknown takes values in the (convex) set 

u) := Qil C R". 

We also assume that a map £ : to — > uniquely determines local equilibria, as 
defined by the relations 

(1.4) QS(u)=u, R(£(u)) = 0, ueuo. 
This suggests to introduce the equilibrium manifold 

M := {U = £{u)}. 

The dimension of the null space of the Nx N matrix B := DR with the equilibrium 
submanifold is assumed to be "maximal" in the sense that 

(1.5) dim (ker(B(£(u)))) = n, 

(1.6) ker (B(£(u))) n lm(B(£(u))) = {0}. 

Finally, since we are interested in the late-time behavior of solutions it is necessary 
to impose that 

QF(£{u)) = c, u £ u), 
for some constant vector c S K™. Indeed, in view of (|1.3p . we can formally write 

ed t u + d x QF{£(u)) ->■ 0, 
so that QF(£(u)) must be a constant, normalized so that 

(1.7) QF = onM 



1 After this paper was completed, P. Marcati pointed out to the authors the relevant references 
|15l 1 16] in which several convergence results are established. 
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Main results. Our objective in this paper is two-fold. First, we determine an 
effective system of equations driving the late-time asymptotic behavior of general 
solutions to (ll.ip in the singular regime e — > 0. Second, we introduce a novel nu- 
merical strategy allowing us to precisely recover the expected asymptotic behavior. 

The paper is organized as follows. Section 2 is devoted to the asymptotic analysis: 
we derive an effective system and, under the hypotheses stated in this introduction 
as well as the assumption of the existence of a mathematical entropy compatible 
with the relaxation source term (in a sense specified below) , we study the structure 
of this system which is found to be of diffusion-type and we show that the associated 
total entropy is non-increasing in time. 

Then, in Section 3 below, we generalize our analysis to a nonlinear version of 
since this is required to encompass certain models arising in the aplications 
when the relaxation has strong nonlinearities. The class of systems under consid- 
eration in this section reads 

(1.8) eSt U + d x F(U) = -W-, 

where the parameter m > 1 introduces an additional scale in the problem. By 
imposing that 

R(£{u) +s U) = e m R(£(u) + MU), U € O, u e w, 

plus certain conditions on the N x N matrix M (see Section 3), we derive again an 
effective system which, now, consists of nonlinear diffusion equations. 

Next, in Section 4, we discuss several specific examples of practical interest and 
show that they fit into our general framework: (1) the Euler equations of compress- 
ible fluids with friction term, (2) the so-called Ml-model of radiative transfer, (3) 
a new model that couples the former two models and, finally, (4) the shallow-water 
system with nonlinear friction. Interestingly, the latter system does require the 
more general formalism (|1.8[) . 

In Section 5, we turn our attention to the discretization of (|1.1[) and (jl.8l) . Recall 
that the "hyperbolic to hyperbolic" relaxation problem was treated numerically first 
by Jin and Xin [22]; our objective is to treat here the "hyperbolic to parabolic" 
relaxation problem. We introduce a new Godunov-type finite volume scheme which 
incorporates a suitable discretization of the source term and allows us to recover 
the expected asymptotic regime. The discrete form of the source term is derived 
by modifying a Riemann solver associated with the homogeneous system. The 
proposed numerical strategy is proven to satisfy a domain invariant principle in 
the sense that all of the computed states belong to the convex set f2. Finally, 
in Section 6, we conclude with several numerical experiments with the proposed 
asymptotic-preserving schemes, and we demonstrate the interest of our method on 
several of the physical models. 

2. Late-time/stiff- relaxation framework 

Derivation of the effective system. Throughout this section, we consider the 
nonlinear system of balance laws (jl.ll) , and we impose the conditions stated in the 
introduction. 

Our objective is to exhibit the system of effective equations satisfied by local 
equilibria u = u{t, x) G lj. In the spirit of Chapman-Enskog expansions of the 
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kinetic theory [8], we consider a formal expansion of solutions U to in the 

form: 

(2.1) U E =£{u)+eU 1 +e 2 U 2 + 

where Ui is referred to as the i th -order corrector. Such an expansion is natural in 
view of the assumptions (|1.4[) and (|1.7p . We consider the system 

(22) ed t W + d x F{U^) = -- £ R{W), 

QU e = u, 

in which we plug the formal expansion (|2.1[) and match together terms of the same 
order of magnitude in e. 

From QU e = u we deduce that Q£{u) = u, and QUi — for all i th -order 
correctors. Next, for sufficiently regular flux F, we obtain the following expansion: 

F(U S ) =F(£(u)) + eA(£{u)) U x + £ l D lF(£(u)) (C/ x , U{) 
+ e 2 A(£(u))U 2 +0(e 3 ). 

Analogously, it is useful to write down the formal expansion of the relaxation term. 
Taking (|1.4p into account, we obtain 

(2.4) -R{U e ) = B{£{u)) Ut + ^-DljR{£{u)).{U u U t ) + sB{£{u)) U 2 + 0{e 2 ). 

£ I 

We then plug (f23|) and ([2^4]) into f[272]) and obtain 
ed t £(u) + d x F{£(u)) + ed x A{£(u)) U x 

= —B{£(u)) U\ — ^DljR(£(u)).(Ui, Ui) — sB(£(u)) U 2 + 0(e 2 ). 



Let us consider (|2.5p and expand in e. The zeroth-order terms give 

d x F(£(u)) = -B(£(u))U 1 . 

At this stage, we can solve for U\ by recalling that QU\ — 0. Indeed, since we have 
imposed the properties (|1.5|) and (|1.6|) on the null space of B(£(u)) then, for all 
fixed J £ Mr, the system 

(2 - 6) = o, 

admits a unique solution V £ K™ if and only if Q J = 0. (See Lemma [2~4l at the end 
of this section.) Here, we have J = —d x F(£(u)) which satisfies Qd x F(£(u)) = 
by (11- 7|) . As a consequence, for all u S w, we can uniquely determine t/i such that 



B(£(u)).tfi = -&F(f («)), 
^• 7j QC/ 1= 0. 

Next, we consider the first-order terms in (|2.5p . and obtain 

ftf (u) + d x A(£(u)) U x = -^DjjR{£{u)).(U u Ux) - B{£{u)) U 2 . 



Multiplying by Q and using the assumption (|1.2p . we obtain 

QDlrR(£(u)).(Ui,Ui) = 0, QB(£{u)) U 2 = 0. 
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Since Q£{u) = u, we arrive at an effective system for the limit u, that is, 

(2.8) d t u=-d x (QA(£(u))U l ), 
where XJ\ is the unique solution to (|2.7[) . 

The role of a mathematical entropy. Now, we assume the existence of a suffi- 
ciently regular, mathematical entropy $ : f2 — > R so that the matrix Dfj$>(U)A(U) 
is symmetric for all U in Q and there exists an entropy-flux map \& : — > R such 
that 

(2.9) Du${U)A{U) = DuV(U), Uefl. 
Hence, all smooth solutions to (|2.2[) satisfy 



(2.10) edt${U £ ) + d x ^{U £ ) = --D U §(U £ )R{U £ ). 

£ 

As usual, we impose <£> to be convex by requiring the N x N matrix Df,Q(U) to be 
positive. In addition, we assume that the entropy is compatible with the relaxation 
in the sense that, for some map v : Ai — > R™, 

(2.11) Du<f>(u)R(u) > o, uen, 1 

(2.12) Du<S> = vQ onM. 

We now analyze the nature of the limiting system (|2 .8[) . 



Theorem 2.1. Consider the nonlinear system of balance laws (jl.ll) . under the 
assumptions (|1.2[) - (11.71) and (|2.1ip - (12.12[) . Then, the associated limiting system 
takes the form 



(2.13) d t u = d x (^SC- 1 (u)S T (d x D u <S>(£(u))) T ^) , 
where 

(2.14) S := QA(£(u)) 
and 

(2.15) C(u):= Dl<S>{E{u))B{£{u)), 

and, for all b £ W N with Qb = 0, £ _1 (m).6 denotes the unique solution to the system 

C{u).V = b, QV = 0. 

In addition, this system is dissipative with respect to the entropy <!> in the sense that 
the following positivity condition holds for all u in uj: 

(d x D u $(£(u))) S£-\u)S T (d x D u <f>(£{u))) T > 0. 

Proof. We follow the strategy in [TT] which we adapt to our problem. We first 
establish (j2.13[) . To simplify the notation, we set 

V{u) = -QA{£{u))U u 

to restate (|2.8p as follows: 

d t u = d x V{u). 

Here, the vector U\ £ R w is the unique solution of (I2.7p . Since D^^(U) is a positive 
N x N matrix for all U € fl, we can rewrite (|2.7I) as follows: 

*(£ (u))i?(£ («)).t^ = -D^(£(u))d x F(£(u)), 

QUx=0. 
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Invoking (|2.15j) . the state vector U\ is defined as the unique solution of 
C(u)U! = -D 2 l} ^{£{u))d x F{£{u)), 
QUx = 0. 

As a consequence, by definition of C~ l {u) we have 

U x = -£- 1 ( U )^$(£( M ))a x F(f (u)). 

From pmjl . we find 

V{u) = SC-\u)I%*{E(u))d m F{e{u)). 
To obtain ()2 . 13[) . we then establish 

(2.17) Dm£(u))d x F(£(u)) = S T v, 
where we define v € K n as follows: 

(2.18) v= (d x D u $(£(u))f . 

To this end we differentiate (|2.12j) (which involves the map v) and we obtain 

(D u £(u)) T Dl<5>{£{u)) =D 2 u <f>(E(u))Q. 
By transposition, we thus get 

(2.19) Dl^{£{u))D u £{u) = Q T 'D^$(f («)). 
Now, we have 

B% $>{£{u))d x F{£ (u)) = ^$(f («))D u f («)fitu 

Since the matrix Dfj<&(£(u))A(£ (u)) is symmetric, we can write 

Dl®{£{u))d x F(£(u)) = {A(£{u))) T Dl$>{£(u))D u £(u)d x u. 
We use (|2.19j) and obtain 

Dl<5>{£{u))d x F{£{u)) = (A(£(u))f Q T D 2 u ^{£{u))d x u 
= S T v, 

and the identity (|2.13p follows. 

The proof will be completed as soon as we establish 

v T SC- 1 (u)S T v > 0, 
which is equivalent to v T T>(u) > where T>(u) — —SUx. We thus have 

v T V{u) = -v T SUi. 
But from (|2~16l) and ([237]) . we deduce 

v T V(u) = (Z#$(£(u))B(£(u))l7 1 ) r [7 1 . 

As a consequence, the expected inequality v T T^(u) > will hold as soon as we prove 
that the matrix Z?^$(£ (u))DjjR(£{u)) is non-negative. 

Recalling the entropy assumption (|2.11[) and the equilibrium property R{£{u)) = 
0, we have 

Du$(U)R(U)>0, Uen, 
(D u $(U)R(U))\ u = £ (u) = 0, u€oj. 
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As a consequence, the matrix Dfj (Du§{U)R(U)) \u—S(u) is non-negative. Next, a 
calculation using the chain rule and the fact that R vanishes on the equilibrium 
submanifold, that is, 

R{£(u)) = 0, 

leads us easily to 

Dfj (Du$(U)R(U)) \ U=£W = D*$(£(u))B(£(u)) + (Dl<$>(£{u))B(£{u))) T , 

in which no third-order derivative terms arise since R precisely vanishes on the 
equilibrium manifold. From the above identity, it then follows that 

(2.21) ((Dfj<l>(£(u))B{£(u)))U) T U > 0, U G ft. 

We thus obtain the expected inequality v Tr D(u) > and the proof is completed. □ 

Monotonicity of the entropy. We then study the asymptotic behavior of the 
entropy inequality. In order to exhibit the entropy law satisfied by the equilibrium 
solution £(u), we need the following technical result. 

Lemma 2.2. Under the assumptions of Theorem \2.1l the entropy-flux map re- 
stricted to the equilibrium ^(£(u)) remains constant for all u inu. 

Proof. We consider the map u i-> ^(£(u)) and, after differentiation, obtain 

£>„*(£(«)) = Du^{£{u))D u £{u). 

By definition of VP given by ()2.9j) . we have 

D u *(£(«)) = Du®{£{u))A(£{u))D u £{u). 

Then, the assumption (12.121) made on $ yields the following relation: 

D u V(£(u)) = D u $(£(u))QA(£(u))D u £(u), 

(2 ' 22) = D u ${£(u))D u QF(£(u)). 

Since we have QF(£(u)) = over ui, then D u QF(£(u)) =0. As a consequence, 
D u ^(£ (u)) = for all uinw and the proof is completed. □ 

Equipped with this result, we can exhibit the asymptotic equation satisfied by 
the equilibrium entropy Arguing the formal asymptotic expansion (|2.1[) 

satisfied by U: 

U £ =£{u) +eUi + ..., 

where U\ is defined as the unique solution to (|2 . T[) , we consider the formal expansion 
of each term in (I2.10[) . First, since the entropy and entropy- flux are regular maps, 
we have 

(2.23) $(U £ ) = $(£(«)) + eDu$(£(u)) XJ X + 0(e 2 ), 
and 

(2.24) V(U £ ) = *(£(u)) + eDu*(£(u)) Ui + 0{e 2 ). 
But, by applying Lemma 12.21 we deduce the following relation: 

(2.25) d x V(U £ ) = ed x Du^(£{u)) Ur + 0(e 2 ). 
Similarly, concerning the entropy relaxation source term, we easily have: 

(2.26) Du^{U £ )R{U e ) = e 2 D 2 I $(£(u))D u R(£{u))U 1 + 0(e 3 ). 
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Now, we plug the expressions (|2.23j) . f)2.25[> and (|2.26[) into (|2.10j) . and consider the 
first-order terms only: 

(2.27) dMS («)) = -8 X (Du*(£(u)) Ux) - U? (D^(£(u))B(£(u))) U u 

where, once again, U± is the unique solution to (|2.7p . From this entropy evolution 
law, we can state an entropy decreasing principle. Indeed, we have established that 
the matrix D\j<$>(£ {u))DuR{£ (u)) is positive (cf. (|2.21|) and, as a consequence, 

u T (dI$(£(u))b{£ («))) £/ > o, c/en. 

We have thus proven the following statement. 

Proposition 2.3. Under the assumptions of Theorem \2.1l the entropy is non- 
increasing in the following sense: 

dM£(u)) < -d x {Du*(£(u)) Ux) . 

To conclude this presentation of the asymptotic system of diffusion equations 
satisfied by stiff relaxation term for late times, let us emphasize the role played 
by the entropy. As recognized by Chen, Levermore and Liu [IT], the existence of 
a convex mathematical entropy provides an important structure to investigate the 
asymptotic regime satisfied by the model. However, the main discrepancy with 
[TTj lies in the nature of the singular limit system. Indeed, in [IT], the singular 
limit system turns out to be an hyperbolic system supplemented by an e first-order 
diffusive term. Here, the obtained asymptotic system defines a system of diffusion 
equations. Even if the limiting solution is smooth (due to the diffusive nature of 
the limiting system), the mathematical entropy is essential in order to establish the 
stability of the asymptotic regime. 

A technical lemma. We end this section with a technical lemma that was useful 
in the above derivation of the asymptotic system. 

Lemma 2.4. Let Abe a real NxN matrix such that dim(ker(A)) = n and kei(A)<~) 
Im{A) — {0}. Let also Q be a real nx N matrix such that rank{Q) = n and QA = 
Then, for any b £ M. N , the linear system 

Ax = b, 

(2 ' 28) Qx = 0, 

admits a unique solution if and only if Qb = 0. 

Proof. If the system (|2.28j) has a solution then left-multiplying Ax = b by Q leads 
to QAx = Qb. Since QA = then Qb = 0. 

Now we suppose that Qb = 0. Since QA — then the columns of Q T are elements 
of the left null-space of A. Furthermore, dim(ker(A T )) = n = rank(Q) therefore 
the columns of Q T are a basis of the left null-space of A. This implies that if 
z e ker(yl T ) then there exists yo € K" such that z = Q T yo- 

Let us consider z £ R. such that Qz = 0. Using the fundamental theorem of 
linear algebra we have z = zq + Zx where zq £ ker(A T ) and zx £ lui(A). 
We then characterize zq and z\ using their respective definitions: there exists yx £ 
M. N such that zx = Ayx, and there exists yo £ M. n such that zq — Q T yo- With these 
characterizations we have: 

= Qz = Qz + Qzx 
= QQ T yo + QAyi = QQ T ya + 0, 
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since QA = 0. Now, rank(Q) = n implies that QQ T is a symmetric positive-definite 
matrix and therefore that yo — 0. Therefore z £ Im(A) if and only if Qz = 0. 

By considering first the existence issue, the above property on b allows us to say 
that b £ lm(A) and hence there exists x £ M. N such that Ax = b. If x is such a 
solution and since kei(A) © Im(A) = l w , we have x = xq + x\ where xq £ ker(A) 
and xi £ lm(A). With these notation xi is a solution to f|2.28j) . 

Considering next the uniqueness issue and suppose that x± and x 2 are two solu- 
tions to (|2.28[) . Then (x\ — x 2 ) is the solution to 

A(xi - x 2 ) = 0, 
Q(xi - x 2 ) = 0, 

and therefore {x\ — x 2 ) £ ker(A) n Im(A) = {0}. Hence, the solution to (|2.28l) is 
unique. □ 



3. Nonlinear diffusive regime 

Derivation of a nonlinear asymptotic system. Some physical models involve 
several relaxation time-scales. More precisely, we suppose now that the ratio of the 
relaxation time and the late time under consideration is no longer constant. The 
extended model thus reads 

(3.1) ed t U + d x F{U) = — -R(U), t > 0, x £ E, 

with U £ fl C R . Here, m > 1 denotes an integer. The case m = 1 has been 
discussed in the previous section, and we now assume m > 1. 

Several of the assumptions made for m = 1 are kept here. Precisey, we assume 
the existence of an n x N matrix Q with rank n < N satisfying (jl.2|) . as well 
as the existence of a map £ : uj — > fl satisfying (1 1 .4[) . We also impose that the 
flux F satisfies (|1.7p . Concerning the nonlinear relaxation map R, an additional 
assumption must be imposed: there exists a N x N matrix, denoted M(e), such 
that 

(3.2) R(£(u)+eU) = e m R(£{u) + M(e)U), U £ Q, uElo. 

The matrix M(e) is assumed to be sufficiently smooth in e £ [0, 1]. 

The assumption initially made on the kernel of the matrix B(£(u)) is irrelevant 
if to > 1. Indeed, this kernel assumption was imposed to ensure the existence and 
uniqueness of the solution to (|2.6|) . We are going to see that this linear system is 
no longer relevant and must be replaced by a nonlinear problem. In this sense, the 
proposed extension is called nonlinear since the diffusive asymptotic regime will 
involve a nonlinear differential operator. 

First, to derive the effective system of equations satisfied by the local equilibrium 
u £ uj, we introduce again a Chapman-Enskog-type expansion: 

U e = £{u) +eU 1 +e 2 U 2 + - 

We plug this expansion into 1|3.1|) and match terms of the same order in e. 

Enforcing QU e = u, the condition (jl.4| on the local equilibrium implies QUi = 
for each corrector term. Note that the expansion (|2.3p for the flux remains valid, 



10 C. BERTHON, P.G. LEFLOCH, AND R. TURPAULT 

and we only have to evaluate the expansion of the relaxation term by recalling (|3.2j) . 
Indeed, this assumption gives 

R(U e ) = e m R(£{u) + M(e)Ui + eM{e)U 2 + ...), 

and so 

j^R(W) = R(£(u) + M(0)Ux) + sB(M(0)Ux).(M(0)U 2 ) 

+ eD e (M(e)U 1 )\ £=Q .R(£(u) + M(0)Ux) +0(e 2 ). 
Setting and $£5$ into ([57T]) . we obtain 

e&£(u) + d^F^u)) + eA(£(u)) f/i 
(3.4) = («) + Af(O)Ui) - eB(M(0)Ux).(M(0)U 2 ) 

- £D £ (M(e)L r i)| e=0 .i?(f (u) + M(0)£/i) + C(e 2 ). 
Considering the zeroth-order terms, we get 

d x F(£(u)) = -R(£(u) + M(0)Ux), 

which turns out to be a nonlinear system of equations with unknown U\ , supple- 
mented by the condition QU\ — 0. 

At this level, we see that the assumption on the kernel of B(£(u)) is no longer 
relevant (or sufficient) in order to solve (|2.6[) . since we now have to consider 

R(£(u) + M(0)tfi) = -d x F(£{u)), 
(3 - 5) QUl = 0. 

In view of the strong nonlinearities involved in this equation, we tacitly assume the 
existence and uniqueness of the solution, denoted Ui, of (|3.5|) . In the applications, 
this property will be checked directly. 

Then, we match first-order terms issuing from p.4[) to get 

d t £(u) + d x A{£ {u)).U x 

= -fl(M(0)Ui).(Af(0)^ 2 ) - D e {M{e)U x )\ e=0 .R{£{u)+M{Q)Ux). 

Since Q£(u) = u for all u£ lj, and since QR(U) = and QB(U) = for all U e O, 
by multiplying the above relation by Q, we obtain 

(3.6) d t u = -d s {QA{£{u))U x ) . 

Once again, this equilibrium equation involves a nonlinear differential operator 
in the right-hand side since U\ is a nonlinear map applied to first-order space 
derivatives. 

Similarly to the "linear" case governed by (|2.8p . we want interpret (|3.6[) as a 
system of diffusion equations. We thus assume the existence of a convex entropy 
$ : 51 — ¥ M which satisfies all the compatibility conditions imposed in the previous 
section. The matrix Dfj^{U)DuF(U) to be symmetric for all U € f2 and, in addi- 
tion, we assume that the compatibility conditions (12. lip and (I2.12p hold. Smooth 
solutions to (|3.ip satisfy the additional balance law 

(3.7) edMU") + d x V(U E ) - -^Du^{U")R{U e ). 
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Equipped with this convex entropy, we observe that U\ is, equivalently, a solution 
to the nonlinear algebraic system 

D^(£(u))R(£(u) + M(0)tfi) = -Dl${E{u))d x F(£{u)), 
QU l = 0, 

for any convex entropy compatible with the relaxation term. 

Next, consider the matrix S defined by (12. 14)) and the vector v given by (|2.18p . 
In view of (|2.20p , the above system is equivalent to 

Dl^{£{u))R{£{u) + M(0)C/ x ) = ~S T v, 
QUi = Q. 

With some abuse of notation and for the sake of clarity, we set 

M U {U{) = Dl§(£(u))R{£{u) + M(0)I7k) 
and introduce the notation 

(3.9) U 1 =N- 1 {-S T v). 
Hence, the equilibrium equation ()3.6j) reads as follows: 

(3.10) d t u = d x (-SAf- 1 (-S T v)). 

Once again, we note the crucial role played by the convex entropies. Indeed, we 
now exhibit the limiting system of equations satisfied by the equilibrium entropy 
$(£(u)). We skip here the details of the computation which are similar to the linear 
case. Lemma l2~2l still holds, as well as the entropy expansion (|2.23j) and (|2.25j) . In 
fact, only the entropy relaxation source term expansion changes. We now obtain 

Du<S>{U E )R(U e ) = e m+1 U?Dtj$(£(u))R(£(u) + M(0)f7i). 

Plugging these expansions into (13. 7p and considering first-order terms, we find 

(3.11) dM£(u)) =-d x (Du^(£ (u))Ui) - UlDl^{£{u))R{£(u) + M(0)Ux). 

To conclude this section, we show that the asymptotic system of equations f|3.10[) is 
of diffusive-type, and that an associated mathematical entropy is non-decreasing. 

Lemma 3.1. Let U\ be given by \3.5)) . Assume the existence of a non-negative 
map c(u) > such that 

(3.12) R(£(u) + M(0)U 1 ) = c(u)U 1 , u e u. 

Then the limiting equation \3. 6)) is nonlinearly dissipative with respect to the entropy 
in the following sense: 

(3.13) -v T SAf- 1 (-S T v)>0, ueuj, 

where S is given by ^2.14\ l and v by \2.18\ l. 

Moreover, the entropy is decreasing as follows: 

(3.14) dM£(u)) < -d x (Du9(£(u)).Ui) . 
Proof. Arguing (|3.8[) and the definition (I3.9p . we have 

-v T SJV-\~S T v) = UfD^(£ (u))R(£ (it) + Af(0)C/i). 
By recalling p,12p . we obtain 

-v T SAf- 1 (-S T v) = c{u)Ul Dm£{u))tJx. 
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The inequality (|3.13|) follows from the convexity of the entropy. Recalling (|3.1ip 
we obtain ()3.14p . and the proof is completed. □ 

4. Physical examples 

Euler equations with friction term. Many models involving distinct physical 
scales enter the framework proposed in the present paper and, specifically, we will 
now illustrate the interest of our framework with four examples. We begin with 
the isentropic Euler equations supplemented with a friction term (cf. [27j [33] for 
further details). The asymptotics for this model has been already considered in the 
literarure and, more recently, relevant numerical techniques have been proposed 
[9]. (See also [26l E].) Importantly, this model satisfies all of the conditions 
required in Section 3, above. 

The Euler model with friction reads 

ed t p + d x pv = 0, 

ed t pv + d x (pv 2 +p(p)) = --pv, 

£ 

where p > denotes the density and v £ R the velocity of a compressible fluid. 
The pressure function p : R + — > R + is assumed to be sufficiently regular and 
satisfy p'(p) > 0, so that the first-order homogeneous system associated with (14. ip 
is strictly hyperbolic. 

In view of (jl.ip , we should set 

(4 - 2) u= (!»)> f(u)= (^+p(p))> R{u) ={pv 

which corresponds to the matrix 

Q = (i o), 

and scalar local equilibria u = p, with 

£{u) = ( J 

As required, we also have QF(£(u)) = 0, and it is easily checked that all our 
assumptions of the previous sections hold. 

Considering now the asymptotic diffusive regime in the limit e — > 0, we note that 
the equilibrium solution must satisfy (12.81) . i.e. 

d tP = -d x (QA(£(u))U 1 ), 

where 

and U\ is the unique solution to (|2.7I) . Since 

=(S d X F(s( U) )=( 9x ° pip) 

the diffusive regime associated with this Euler model with friction is described by 
the equation 

(4.3) d t p = dl{p{p)). 

Based on Theorem 12. 11 we observe that the diffusive nature of (|4.3p follows from 
the existence of a convex entropy which is compatible with the relaxation source 
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term in the sense (12.11[) - (|2.12[) . Indeed, by introducing the internal energy e(p) > 
defined by 

P(P) 



e'{p) = 



P 2 



we see that smooth solutions to (I4.1[) satisfy 

£dt ( p y + + dx (p~2 + P e (p) + p(p^j v = ~\ pv2 - 

2 

The function $(f ) = p\ + pe(p) is a convex entropy satisfying the compatibility 
conditions with the relaxation terms. 

The Ml model for radiative transfer. The second example of interest relies 
on a more complex physical set-up, relevant in radiative transfer and referred to as 
the Ml-model [IH [29]. (See also [6j ED].) This model reads 

ed t e + d x f = i(r 4 -e), 

£ 

(4-4) ed t f + d x (x{f/e)e)=--J, 

ed t r = ~(e - r 4 ), 
e 

where e > is the radiative energy and / the radiative flux, restricted by the "flux 
limitation" condition 

'- <i. 

e 

and r > denotes the temperature. The function \ '■ [~ 1> 1] — ^ K + stands for the 
Eddington factor defined by 

5 + 2^4^3^ 

Once again, we rely on the notation introduced in the previous sections (cf. sys- 
tem Ql.ip ) and we write 

F(U)= | X(i)e | , R(U)= | / 
The local equilibria are described by the map 

£{u) = 

where u = r + r 4 is now a scalar. We set Q = (1 1) and, in agreement with (jl.7l) . 
we have QF(£{u)) = 0. 

The asymptotic regime is governed by (|2.8j) and to exhibit its explicit formulation 
we need the expression of A(£(u)) and U±. A straightforward calculation gives 

A{£(u)) = 







1 






1 




x(o) 


x'(o) 


: ■ 

















\l 
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while Ui is the solution to (12.71) which here reads 




We thus obtain 



and, according with (|2.8|l . the asymptotic regime is governed by the following dif- 
fusion equation: 



A coupled Euler/Afl model. We propose here an example of system that de- 
generates into a system of diffusion equations of dimension n > 1. To do so, we 
couple the Euler model (|4.1[) with the Ml model (|4.4|) as follows: 

ed t p + d x pv = 0, 



(4.5) 



edtpv + d x (pv + p(p)) = — pv + -f, 

e e 

ed t e + d x f = 0, 
7 



ed t f + d x x 



e=-~L 

e 



in the notation previously introduced. Here, k and a denote positive constants. 
Even though this is a toy model, the pressure has to be sufficiently small in order 
to represent an application of physical interest. We will therefore consider the 
following pressure law: 



p(p) = c vP \ 

In the formalism (jTTTJ) we need to set 



C p < 1, 



U = 



F(U) = 



pv 

e 

The local equilibrium is given by 




pv 2 +p(p) 



\ 



T) > 1. 



R(U) 



V X(|)e / 



/ o 

KpV — of 



V 'f J 



8{u) 



e 



u = QU = 



Q = 



10 
10 



Once again, one has QF(£(u)) — 0. To derive the asymptotic regime, we exhibit 
A{S{u)) and Uf. 

( o 

\ [-dxP(p) ~ \d x e] 



A(S(u)) = 



/ 1 0\ 

p'{p) 

1 

V o o i o/ 



Ui 



V 







-4zd x e 



J 
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The asymptotic diffusive regime of the system (|4.5|) is therefore given by 

dtp - -d 2 x p{p) - ^-d 2 x e = 0, 

(4.6) K 6K 

d t e - —d 2 x e = 0. 

OCT 

Shallow water with strong friction effect. The last suggested example is de- 
voted to the well-known shallow-water model supplemented by a strong friction 
term in a late-time regime where the friction is assumed to dominate over the 
convection. 

This model is given as follows: 

ed t h + d x hv = 0, 

(4.7) K 2(h) 

ed t hv + d x (hv 2 +p(h)) = -^-ghv\hv\, 

where ft. > is the water height, v € R the velocity and p{h) = g^- the pressure law. 
Here, g > is the usual gravity constant, while the friction coefficient k : M + — > R+ 
is a given and positive function. In [28], several examples of friction k are proposed. 
A standard choice is n{h) — ^ where kq > is a given parameter. 

We note that this model enters the framework of the nonlinear extension gov- 
erned by (|3.1[) with m = 2. According to the notation involved in (|3.ip . we have 
set 

U= (hv)-> FiJJ) = ( hv 2 +p(h) ) ' R{JJ) = ( K 2 (h)ghv\hv\ 
Concerning the equilibrium, we have 

' h 


where u — h is a scalar and Q = (1 0). The assumption (|3.2[) easily holds since 
R(S (u) + eU) = s 2 R(£{U) + M(s)U), 

where 

We turn our attention now to the asymptotic regime which is governed by the 
nonlinear diffusion equation (|3.6p . To get its explicit from, we have to exhibit 
U\ = (a P) T the solution to (|3.5p which reads as follows: 

a — 0, 

K 2 (h)g(3\(3\ = -d xP (h). 

We easily find 

_ \fhd x h 
K.{h)yf\d x h\ ' 

and the effective nonlinear diffusion equation is thus 
(4-8) d t h = dj^- -^L=) . 

This is a nonlinear Laplacian equation (for instance, see |23j and references therein) . 
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Lemma 13.11 applies here with the following entropy. By introducing the internal 
energy e(h) — gh/2, we see that smooth solutions to (|4.7|) satisfy 

edt (^y +5y) +d x (h^+ghA v = -^-ghv*\hv\. 

2 7 2 

The entropy $(J7) = h\ + satisfies all the required properties, and the condi- 
tion (|3.12p holds since 



R(£(u) +M(0)J7i) = 




where Ui = (0 fi) T . As a consequence, we obtain R{£{u) + M (0)U\) — c(u)U\ with 

c(u) = gn(h)\/h\d x h\ > 0. 
Hence, the limit equation (|4.8p is of diffusive- type in the sense of Lemma [37X1 

5. Asymptotic-preserving schemes 

Objective. In this section, we consider the numerical approximation of solutions 
to (jl.ll) . Our goal is to derive a class of numerical schemes that restore the relevant 
asymptotic regime, given by (|2.8|) . in the limit e — > 0. One of the main difficulties 
when deriving asymptotic-preserving schemes lies in the independent role played 
by each e and the mesh size. More precisely, the limit discrete diffusion equation 
(as e tends to zero) must be obtained independently of the space mesh-size. 

Such a numerical problem was investigated during the last decade on several 
specific examples. For instance, in [5], the Euler equations with friction term are 
considered. Relating works to radiative transfer and the Ml-model are given in 
O [20] . The reader is also referred to [3] where distinct physical applications are 
proposed. 

In the present work, we propose a generalization of a numerical scheme derived to 
approximate the solutions to the Euler equations with friction and the Ml-model 
[3J. To sketch this suggested numerical procedure, we first consider a standard 
finite volume method to approximate weak solutions to the homogeneous system 
associated with (jl.l[) in which we have omitted e: 

(5.1) d t U + d x F{U) = 0. 

Next, we derive a suitable correction to obtain a finite volume discretization of 
the source term. Hence, the corrected finite volume method gives approximate 
solutions of the following system: 

(5.2) d t U + 3 x F(U) =-<yR(U), 

where 7 > is a fixed parameter. Finally, the asymptotic behavior of the scheme 
is analyzed. We consider a late-time compatible discretization and we fix 7 = -, 
The asymptotic scheme is thus obtained in the limit of e to zero. Modulo a suitable 
correction, it is thus proved to be asymptotic preserving. 

A discretization of (|5.2[) . As a first step, we suggest to consider the well-known 
Godunov-type scheme introduced by Harten, Lax and van Leer [21] with a single 
constant intermediate state, to approximate the weak solutions to (|5.1|) . 

We consider a uniform mesh made of cells [^1-1/2)^1+1/2) where x i+ i/ 2 — Xi + 4p 
for all i € Z with a constant cell size Ax. The time discretization is defined by 
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t n+1 = t n + At where the time increment will be restricted later on by a CFL like 
condition. 

We define the discrete initial data as follows: 

U (x) = — / U(x,0)dx, x € [x i - 1 /2 ) x i+1 / 2 ). 

We seek for a piecewise constant approximation of the exact solution of (|5.1[) at the 
time t n , 

U n (x) = U?, 16 [x^ 1/2 ,x, +1/2 ), 
with UJ 1 E Q, for all i E Z. 

By considering a suitable sequence of approximate Riemann solvers, we can 
evolve this approximation and get a piecewise constant function U n+l {x) which is 
an approximation of the solution to (15 . 1[) at the time t n + At. Following Harten, Lax 
and van Leer |21] , at each cell interface we use the following approximate Riemann 
solver: 

U L , ?<-&, 



(5.3) 



Un(-;U L ,U R ) 



t 

U\ -b<-<b, 
t 

Ur, - t >b, 



where b > is a fixed and sufficiently large constant, and 

(5-4) U* = \{Ul + U r ) - ^(F(U R ) - F(U L )). 

As a consequence, as soon as the following CFL restriction holds: 

(5.5) 



b At < 1 
~Ax ~ 2'' 



we are considering a juxtaposition of non-interacting approximate Riemann solver 
(cf. FigureE]) denoted Ug x (x,t n + t) for t £ [0, Ai). 

-b b -b b 





x i-l/2 Xi+1/2 

Figure 1. HLL scheme: Juxtaposition of approximated Riemann problems. 



The updated approximated solution at time t n+1 is thus defined as follows: 
(5.6) 

Since we have 

fl* , = 

2 V "' ' ~ l+u 2b' 



Ul l+1 = ±- Ul x (x,t n + At)dx. 



(5-7) U? +1/2 = hu? + u: +l ) ±(F(U? +1 ) F(U?)), 
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an easy computation gives the following standard conservation form: 

— (F HhI f - F HL1 T ) 
Ax i+1/2 i- 1 / 2 ' 



(5-8) U- +1 =U-- ^-{F^A - F^ 2 ), 



(5-9) F^f 2 = X -{F{U?) + F(U? +1 )) - "-{U? +l UT 



where we have set 

I(F(C/f) + F(C/f +1 ))-^( 

For simplicity in the presentation, we have adopted here a constant numerical 
cone of dependence (cf. Figure Q]) characterized by a single speed parameter b > 0. 
As prescribed in [2T] (cf. also [2H 132]), each cone of dependence can be variable 
and defined by a pair (b~ +1 / 2 , bf +1 , 2 ) with b~ +1 ^ 2 < bf +1 , 2 . However, for the sake 
of simplicity and without genuine loss of generality, we present our strategy for the 
simpler case bf +1 , 2 = —b^ +1 , 2 = b > 0. 

At this level, a first remark must be done concerning the invariant domain prop- 
erty for the scheme (|5.8j) . It suffices to ensure that U* +1 , 2 belongs to f2 for all i G Z 
to deduce that £/™ +1 in f2 for all i € Z. Indeed, from (|5.6[) . we have 



and, based on the CFL restriction (I5.5[) . the above relation is a convex combination 
of states in Q. Since 51 is a convex set, we deduce that U™ +1 belongs to Q. 

Our main necessary condition for the above argument is that U* +1 , 2 £ O for all 
i in Z. But, once again, U* +1 ^ 9 can be seen as a convex combination, as follows: 

U? + i/a = \ fa + - b F(un) + \ (Vf+i - \f{U? +1 ) 

Since Vt is an open convex set, we can choose b to be large enough so that to enforce 
the condition U* +1 / 2 € Si- 
Next, we modify this approximate Riemann solver and introduce a discretiza- 
tion of the source-term in order to approximate the solutions of (|5.2p . Similar 
modifications were made for specific problems in [5] and p, 6 , while we propose 
here a general approach based on matrix-valued free-parameters. We modify the 
approximate Riemann solver (|5.3[) as follows: 



„ x 
Ul, - < 



(5.10) U n (j;U L ,U R ) = { 



t 

U* L , -b<-<0, 
t 



U* R , < - < 6, 
U R , I > 6, 



where we have set 
(5.H) 



U* L = all* + (I - a)(U L - R(U L )), 
U* R = all* + a)(U R - R(U R )). 
Here, a denotes a N x TV matrix defined as follows: 



M->) il= (1+^(1 + 2.) 
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R(U) = (I + a^RiU). 



and 
(5.13) 

The NxN matrices I and a respectively denote the identity matrix and a parameter 
matrix to be defined. 

A choice for the matrix a will be made later, which will turn out to govern the 
asymptotic regime. At this point, the matrix a is assumed to be such that the 
inverse matrices in (|5.12|) and (|5.13j) are well-defined. 

We adopt the modified approximate Riemann solver (|5.10[) to derive a modified 
Godunov type scheme in the spirit of [3]. At each cell interface 2^+1/2, we set the 
approximate Riemann solver U-r ( x 1 2 ; £/" , U™ + 1 ) to define a juxtaposition of 
modified approximate Riemann solver, denoted by U£ x (x,t n + t) for t G [0, At). 
(See Figure El) 

-b b -b b 





x i~l/2 x i+l/2 

Figure 2. Juxtaposition of modified approximated Riemann problems. 



Thanks to the CFL condition (|5.5p . such a juxtaposition is non-interacting. 
At time t n+1 , the updated approximated solution is given as follows for all i € Z: 



rn+l 



(5.14) UP + = / UAx(x,t n + At)dx. 

We compute this integral form and obtain 

— (U" +1 - [/") + , 4 I 2 ' ii— _> ' , I 2 ' 

-^('-Oi+l^^wm 

where the numerical flux F^y 2 is given by ()5.9|) . The discretized source-term can 
be rewritten in a more relevant form: 

- a i+1/2 )R z+ i/ 2 (Ul l ) = -^a. l+1/2 (a^ 1/2 - I)Ri +1 / 2 (U" ). 

In view of the definition of Oi+1/2 (see (15.12[) ). we deduce that 



A 



(7 - 0^/2)^+1/2(1/?) = ^+1/2 



and, similarly, 
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As a consequence, the scheme (|5.15|) reads 
(5.16) 



This scheme satisfies the following statement. 

Theorem 5.1. Assume that the matrix a defines a spatially continuous map. The 
numerical scheme H5.16\) is consistant with the equation \5. 2\) . 

At time t n , assume U" € fl for all i 6Z and, in addition, that the state vectors 
an d defined by 

= a l+1/2 U? +1/2 + (I - a l+1/2 )(Ul l - R(U?)), 

Ui+l/2 = "i+1/2^+1/2 + i 1 ~ m+l/2)(Ui + l - R{Ui+l)), 

belong to il. Then, the updated state vector t/™ +1 , defined by i5.16}) . belongs the 
set f2 for all i in Z. 

Proof. The consistency property follows from the definition of a i+1 i 2 , given by 
(|5.12jl . Indeed, we easily obtain 

a i+1/2 = r + o(Ax), 

to deduce the expected consistency of the flux and the relaxation source term. Only 
the term 3^(^+1/2 ~ ^i~i/ 2 ) F {Ui) is left over. Recalling (|5.12[) . we have 

1 7 

^(fii+l/2 - a i -l/ 2 ) i? (C / ") = -^7ai+l/2fc+l/2 - °4-l/2m-l/2 F ( U i )> 

to obtain 

^fe+1/2 ~ Qk-i/iMU?) = O(Ax), 

as soon as g_ i+ i/ 2 — 24-1/2 = O(Ax). The expected equation consistency is therefore 
obtained. 

Concerning the robustness of the method, from (|5.14p . we have 

r«+l - h—TT+L , {■, _0h — \ TTn ,h—TT*R 



U < = b A~ X U ^ + y 1 - 2b Ax~) W + b A~ X U ^ 

which is nothing but a convex combination in J7. Then U™ +1 is in Q and the proof 
is completed. □ 



To conclude this derivation, observe that the term ^(Si+1/2 — £k-i/2)F(U™) 
may seem to be a discrepancy in the method. In fact, this term is standard to 
derive asymptotic preserving schemes and it can be found in several works. (See, 
for instance, [U[3J[3].) Our approach allows us recover a scheme proposed earlier 
in [TS]. 
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The linear asymptotic regime. The scheme (|5.16[) is now considered to suggest 
a discretization of our initial model (jl.lj) . The expected scheme is thus easily 
obtained when substituting At by =^ and fixing 7 = ~. The resulting scheme reads 
as follows: 

±. {U n + i _ un + ±{a i+ y 2 F^-^_ y2 F» L j 2 ) 



(5.17) 



Ax. ' ' 



where 

(5-18) S*+i/2= ^+^(-f + £i+i/ 2 ) 

For the sake of simplicity in the forthcoming asymptotic derivation, we propose to 
introduce the N x N matrix o| +1 y 2 defined by 

f Ax N _1 

(5-19) o^ 1/2 = \ eI+—(I + a i+1/2 ] 

so that we have a i+1 / 2 = £a? +1 / 2 . Recalling the definition, the scheme (|5.16j) takes 
the form: 

(5.20) At Aa; , 

We observe that 17™ remains close to the equilibrium state for e small and 

we thus consider the following expansion: 

C/f = £(<)+ 6(1^)? + G(e 2 ), 
which we now plug into (|5.20p . We easily have 

-R(Un=B(£(u2)).(U 1 )- + 0(e), 

£ 

pHLL 1 _ [.ALL 



where 



^i+l/2\£(u)+0(e) ~ F i+\/2\E{v) + C( £ )> 



F"ifMu) = \ (F(£K)) + F(£K +1 ))) \ {£ « +1 ) - £ «)) . 



In addition we have 



f i+ i,2 = ^-{I + <h + x,2T l +0{e). 



Ax 

By considering the first-order terms issuing from (|5.20p . we obtain 

^(£« +1 ) -£(<)) 

+ (( 7 +2i+i/2)~ 1 - F i+i/2l£:(u) - U + 2i-l/2) _1 - F i-l/2|£(«)) 

= ^ ((/ + ^i/a)- 1 - (/ + ^-i/a)" 1 ) 

- ^ ((J + ffif+l/2)" 1 + U + ^-l/2)- 1 ) 
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Next, assume the existence of a n x n squared matrix, denoted by .M j+1/2 , such 
that 

(5-21) Q(/ + ^ +1/2 )- 1 = ^X !+ i/ 2 0. 

Recalling the assumptions (|1.2p . (|1.4p . and (|1.7|) . we write 



HLL I 

-l/2\£(u) ) ; 



where 

Q*£v 2 k«) = (m<)) + m<+i))) - \q - £(<)) 

As a consequence, the asymptotic discrete regime is given by 

(5.22) - <) = J-j (7W i+1/2 « +1 - <) + M-1/2 W-i - <)) . 
We thus have established the following result. 

Theorem 5.2. Consider any N xN matrix g_ i+1 ^ 2 such that the matrices I+o_i + i/ 2 
and (1 + + 24+1/2 are nonsingular for all £ > 0. Assume the existence of a 

nx n matrix such that 115.21]) holds, and introduce the nx n matrix A4(u) 

defined by 

M{u) = QA{E(u))L' 1 {u)Du^{E{u))A(£{u))D u E(u). 

In addition, assume that the matrix Mi+x/i is a discrete form of M(u) at each 
cell interface x i+ i/ 2 . Then, the asymptotic behavior of the scheme \5.2U\) coincides 
with a discrete form for the limit diffusion equation 12.13]) . 

Proof. We directly deduce from (|2.13j) that the diffusive limit equation reads 

d t u = d x (M{u)d x u) . 

Since the asymptotic regime satisfied by the scheme (|5.20p is governed by (|5.22[) . 
the proposed choice of the matrix M i+ i/ 2 leads us to the correct behavior of the 
scheme as e goes to zero. □ 

The nonlinear asymptotic regime. We propose to extend the above numerical 
scheme to consider the nonlinear asymptotic regime governed by the system (13. ip . 
To address such an issue, once again we consider the scheme (|5.16p where At is 
substituted by — and we set 7 = ~. To be consistent, we substitute R(U- 1 ) by 

Adopting such a strategy, the same arguments used to obtain (|5.20|) now give 

i-(t/f +1 - U?) + ^-Kfi/a^va " ^-1/2^/2) 

(5.23) At Aa; 

= ^tern/a - «?-i/ 2 )nE?) ^ ^fern/a + a!_ 1/2 )R(U^). 

where cf i+1 / 2 is defined by (|5.19p . The above scheme exactly coincides with (|5.20p 
as soon as we fix m = 1. 
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Once again, as soon as e tends to zero, the state vector JJf remains in a neigh- 
borhood of the equilibrium state £(u"). As a consequence, we adopt a formal 
expansion given by 

U? =£«)+ s(W + 0(e 2 ). 

Arguing the same calculations as used in the linear expansion case, the asymptotic 
discrete equation is given by (|5.22[) . Hence, we have to choose Ai i+ i/ 2 in order to 
get a discretization of (13.61) . 



Theorem 5.3. Consider an N x N matrix 24+1/2 suc h that the matrices 1+24+1/2 
and (1 + 2^)/ + 2i+i/2 are non-singular for all e > 0. Assume the existence of a 
n x n matrix M. 1+1/2 such that \5.21)) holds. Consider also the n x n matrix M.(u) 
defined by 

(5.24) M(u) = QA(£(u))K~ 1 (u), 

where IZ^ 1 : cj — > Q, defines the unique solution of \3.5\) . Assume that M-i+i/2 * s 
a discrete form of M.(u) at each cell interface x i+ i/ 2 . The asymptotic behavior of 
the scheme h5.23)) defines a discrete form of the nonlinear diffusion equation hS.b}) . 

Proof. By definition of 7?. _1 (m), we have U\ = 7£ _1 (u) the unique solution to (|3.5|) . 
As a consequence, we deduce the following rewriting of (|3.6j) : 

d t u = d x {M{u)d x u) , 

where M.(u) is given by (I5.24[) . The proposed definition of M.1+1/2 concludes the 
proof. □ 

6. Numerical experiments 

Euler equations with friction. To illustrate the interest of the asymptotic pre- 
serving scheme (|5.17p , we apply it now to the Euler equations with a friction term 
(|4.ip . Here, we suggest to fix the matrix parameter 24+1/2 as follows: 

94+1/2 = cr i+l/2^ J 

where <Xj-|-i/2 stands for a scalar parameter to be defined. As a consequence, the 
matrix a i+1 / 2 , defined by (|5.18[) . is now given by 

1 



(6-1) a i+1/2 = a i+1/2 I, a t+1/2 = A r. 

1 + 355( 1 + <7 <+i/a) 

The scheme f|5 . 1 7[) thus reads 



(6.2) 



e / n+l n \ 1 I t^P,HLL rip,HLL 

^{Pi - Pi) + (a i+ i/2if+ 1/2 - Oi_ 1/2 F£_ 1/2 



Vi+1/2 - Vi-1/2 , .„ 

-o.i+i/2 "i-i/aC/w)* , 



(6-3) = - ai+1/2 a J±^-pzlll a i ^/2 (p?W? +p(p?)) 

e 2 

where the numerical flux vector '^i+i/2 ) * s defined by (|5.9|) . 
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First, applying Theorem 15. 11 we observe that the proposed scheme is consistent 
with the system (|4.1|) and preserves the admissible states. To address such an issue, 
in view of Theorem 15. 1[ we have to establish the positivity property (i G Z) 

Pi+l/2 = a t+l/2Pi+l/2 + (1 _ a i+l/2)P™, 
P*+l/2 = a t+l/2Pi+l/2 + i 1 - a i+l/2)Pi+l- 

Since a i+1 / 2 , defined by (|6.1[) . belongs to (0, 1), we have P*^_ X I2 > ® an d Pi+i/2 > 
as soon as p* +1 / 2 > which is satisfied for sufficiently large values of b (cf. relation 

(El). 

Now, we study for the asymptotic behavior of the scheme (16.2[) - (|6.3[) . Put in 
other words, we have to fix the free parameter cr i+1 / 2 to recover the expected 
asymptotic regime in the limit of e to zero. This required asymptotic behavior 
must be governed by a discrete form of (|4.3[) . First, observe that 

a i+l/2 

lim ati-L.\/2 = 0, lim = 

e^O ' e-K) £ 

From the momentum (pv)™ +1 governing equation (|6.2p . we easily deduce the fol- 
lowing momentum behavior in the limit of s to zero: 

( P v)2 = o, i g z. 

The density approximation thus admits the following asymptotic regime: 

1 (n n+1 n n \4- 2b ( 1 F P,HLL, 1 ppMLL, \_ n 

Ai^i - Pi ) + ^2 ^ 1 + C7 . +1/ /i + i/ 2 U=o - i + ^/i-i/a ^=°J - °- 
But we have 

*i-i/ 2 \pv=o - ~2 lP*+i 



1 , „_ui & 2 / I 



and therefore 
We choose 



1+1 1 - 1 if n™ ^ rt n 

i ) 11 Pi+1 r Pi ) 



(6-4) CT l+ l/2 



P(P?+l) -P(P?) 
b 2 



1, otherwise, 



P'(p? 

and arrive at the following discretization of the diffusion equation 

^(p- i+1 - p?) = (p(p?+i) - MpD +Kpti)) 



To conclude, observe that Theorem 15.21 applies if the matrix Mi + \/ 2 is defined by 
■M-i+i/2 = m i+i/2^n where /„ denotes the n x n identity matrix and 



Kp"+i)-Kp?) 



ifp?+l^p", 



rrn+i/2 = \ Pi+i ~ Pi 

p'(Pi), otherwise. 

In [2112], the scheme (|6.2[) - (|6.3p was derived by a completely different approach. 
Similarly, in [3J [21 E] , the application of our general asymptotic preserving scheme 
(15.171) is found in the framework of the radiative transfer (I4.4[) . 
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Figure 3. Reference (full line) and proposed scheme (dashed line) 
solution comparison. 



'Entropy - 600 cells' 
'Entropy - 300 cells' 



Figure 4. Decrease of the entropy. 



As an illustration, the scheme (|6.2[1 - (16.31) supplemented by the asymptotic pre- 
serving correction (|6 .4[) is used to approximate the solution when the initial data 
is given by 

rp J(2,0) T , x G [1.2, 1.8], 
I (1,0) , otherwise. 

Furthermore, we choose the simple pressure law p(p) — p 2 and e = 10 -3 . In 
Figure [3l the solution on a 300 cells grid (Ax = 10~ 2 ) is compared at time t = 



2(5 
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2.10 4 with a numerical approximation of (|4.3[) computed with 600 cells. We note 
a fairly good agreement between the two approximate results. In agreement with 
Lemma 12.31 Figure U shows that entropy is decreasing in time. For two choices of 

space grids, we plot Y1\P\ + P e (p)j. versus time. 



Coupled Euler/Afl equations. We now propose to approximate the solution 
of the system (|4.5[) by adopting the asymptotic-preserving scheme (|5.17p with the 
following matrix parameter cr.j +1 / 2 : 



y+i/2 



(°-i.i+i/2 ~Sa,i+i/2 °\ 




V o 



23,i+l/2 











0/ 



where g_j f+j/g & re parameters to be defined later, in order to reach the required 
asymptotic regime (|4.6p . With such a definition a| +1 ; 2 , defined by ()5.19j) . becomes: 



(2b k e 



\ 











2b k e 



2b k e r yAxa_ 



2,k 



2bkS + jAx 




o 

2b k e 







"3, He 




o \ 




2b k e 
2b k e + ^AxJ 



where 8j tk = 2b k e + jAx(l + g_j k ) and 7 = 
associated with this scheme is 



max(K,<r) 



The asymptotic regime 



„n+l n 
Pi = Pi 



At b 2 



Pi+i ~ Pi 



n _ n 
Pi Pi-1 



Ax j Ax yi + £ 1)i+1/2 l+£i,i-i/ 2 

£2,i+l/2( e i+l — e ?) 



<L2,i-l/2( e ? e ?-l 



(1 + ct m+1/2 )(1 + £3^+1/2) (1 +£ M _i/ 2 )(l +£3,1-1/2) 7 ' 



n+1 r. 

e„- = e„ 



'i+1 



At fe 2 



Az 7A2: \1 +£ 3 , i+ i/ 2 



1 + <T 



■3,2-1/2 



To recover a discrete form of (|4.6p , we propose to set 



£1,1+1/2 



12,i+l/2 



-3,i+l/2 



« gj+1 - Pi 

7 P^+i - P? 

3(7 



1. 



7 
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In conclusion, setting Ap™ +1 , 2 



P'i 



+ 1 



J i+l 



24+1/2 



V 



T A pr+i/2 
o 






Pi 

■ 1 



we obtain 








T A P?+l/2 



3a_ 1 

7 










This scheme is now applied to the following numerical experiment. We choose 
an initial data given by 



0.2, 



0, f = 0, 



1.5, x e [0.45,0.55], 
1, otherwise. 



1, e 



10' 



and r) = 2. 



The parameters of the model are n = 2, a 
The numerical solution is computed with Aa: = 10 -2 and compared to Figure[5]with 
a reference solution obtained by solving (|4.6|) . Once again, the solution is perfectly 
captured even on a coarse grid. We emphasize that this case is very challenging 
because of the very specific form of diffusion involved on the density. 

7. Concluding remarks 

We have presented a general framerwork to investigate the late-time/stiff-rela- 
xation limit of a large class of hyperbolic systems. This framework was shown to 
cover the examples of interest arising in the modeling of complex fluid flows when 
several time-scales are involved. A new class of schemes was proposed for their 
approximation, and we demonstrated that the proposed modification was crucial 
in order to ensure the correct asymptotic behavior for late times. It would be 
interesting to make comparisons between the numerical results obtained here and 
concrete experimental results, especially on radiative transfer problems. In future 
work, we plan to generalize our continuous and discrete frameworks to problems 
with several space dimensions, and to numerically investigate the robustness and 
accuracy of such multi-dimensional finite volume discretizations in realistic physical 
situations. 
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